An important problem in the field of electromagnetics is the analysis and design of radiating and scattering objects. In one formulation, the electric and magnetic fields scattered/radiated from a material object are given by:
                    E        _            s        =                  ∮        s            ⁢                        [                                                    jωμ                (                                                      n                    ⋒                                    ×                                      H                    _                                                  )                            ⁢              ψ                        +                                          (                                                      n                    ⋒                                    ×                                      E                    _                                                  )                            ×                              ∇                ψ                                      +                                          (                                                      n                    ⋒                                    ·                                      E                    _                                                  )                            ⁢                              ∇                ψ                                              ]                ⁢                  ⅆ          S                                        H        _            s        =                  ∮        s            ⁢                        [                                                    jωɛ                (                                                      n                    ⋒                                    ×                                      E                    _                                                  )                            ⁢              ψ                        -                                          (                                                      n                    ⋒                                    ×                                      H                    _                                                  )                            ×                              ∇                ψ                                      +                                          (                                                      n                    ⋒                                    ·                                      H                    _                                                  )                            ⁢                              ∇                ψ                                              ]                ⁢                  ⅆ          S                    where Es and Hs are the scattered/radiated electric and magnetic fields. The vectors, E and H, are the total electric and magnetic fields. The tangential and normal components of the total fields are interpreted as currents and charges: J={circumflex over (n)}× H electric current  M={circumflex over (n)}×Ē magnetic current ρ=ε{circumflex over (n)}·Ē electric charge ρ*=μ{circumflex over (n)}· H magnetic charge The function ψ is a Greens function that may be interpreted as an outward-traveling spherical wavelet:
  ψ  =            ⅇ              j        ⁢                                  ⁢        kR                    4      ⁢      π      ⁢                          ⁢      R      where R is the distance from a source point to an observation point. The scattered fields at an observation point are thus viewed as the summation of spherical wavelets emanating from electric and magnetic current sources and charges on the boundary of the scattering object.
By application of boundary conditions, the above-written integral equations can always be transformed, by means known in the art, into a system matrix equation of the form:V=ZJwhere V is a known vector or matrix excitation function, Z is a calculable system matrix, and J is an unknown source vector or matrix to be determined which can represent electric and/or magnetic currents. A matrix element of Z, Zij, quantifies the electromagnetic interaction between an ith source point and a jth observation point. This interaction is characteristic of the object, independent of the value of the source function at the source point or the field at the observation point. Generally, a source value, Ji, at a source point i, depends on the excitation values, Vj, at many observation points, j. Stated conversely, the excitation Vj at an observation point, j, depends on the source values Ji at many source points, i. In the usual case, the values Vj and Zij can be computed and the system matrix equation is solved for each Ji. Once the Ji are determined, the unknown fields scattered or radiated from the object may be computed directly from the integral equations given above.
For purposes of clarity we will show how the solution is obtained for the case of scattering from a perfectly conducting object. For a perfect electrical conducting object, magnetic currents and magnetic charges are zero. The boundary conditions are:{circumflex over (n)}×Ē={circumflex over (n)}×(Ēs+Ēi)=0{circumflex over (n)}× H={circumflex over (n)}×( Hs+ Hi)= Jwhere the tangential components of the incident fields, n×Ei and n×Hi, on the surface of the body are specified. Using these results along with the continuity equation for charge conservation, the electric field integral equation above reduces to:
            n      ^        ×                  E        _            i        =            n      ^        ×          ∫                        [                                    jωμ              ⁢                                                          ⁢                              J                _                            ⁢              ψ                        +                                          j                ωɛ                            ⁢                              ∇                                  ·                                      J                    _                                                              ⁢                              ∇                ψ                                              ]                ⁢                  ⅆ          S                    Thus, the problem is to determine the electrical currents, J, on the body, given the tangential component of the incident electric field, n×Ei, on the body. In this formulation, n×Ei is the excitation function and J is the source function.
To solve for the unknown electrical currents, the integral equation is discretized into a matrix equation. To accomplish this, the current, J, is expanded in a set of basis functions. The chosen basis functions depend on the geometry of the object and how that geometry is modeled. Thus, for example, the object may be two or three dimensional and modeled with wire mesh or with surface patches that may be rectangular or triangular. (A typical triangle basis function for surfaces is the Rao-Wilton-Glisson (RWG) basis function. This comprises two surface triangles joined by a common edge. The vector part is from the center of one triangle to the center of the second triangle and is the vector direction.). Then, the integration and differential operators can be computed on the source coordinates at each surface patch. To complete the discretization process, both sides of the integral equation are multiplied by a vector of weighting functions and integrated over the observation coordinates. These techniques outlined above are known by persons of skill in the art and the method is known generally as the Method of Moments (MOM).
As noted, the discretization of an integral equation produces a system matrix equation of the form:V=ZJwhere V is a known excitation function, Z is the system interaction matrix which can be calculated, and J is a source function to be determined. Note that the matrix Z is characteristic of the object and does not depend on the source or excitation. Once J is determined, it can be substituted into the original integral equation to solve for the scattered/radiated field. In particular, one may compute a far field disturbance arising from the sources, J. If the matrix Z were small, then J could be determined by:J=Z−1VFor objects of electrically moderate or large size, the matrix Z is very large. The time it takes to compute the elements of Z and to compute the inverse of Z is very long. For a problem of practical size, inverting Z could take many days even using an exceptionally fast computer. Moreover, due to limits on machine precision, the solution obtained computing the inverse is insufficiently accurate. And for many computers and problem sizes, Z is too large to be inverted at all. For these reasons, solving the equation by inverting the system interaction matrix is probably never done for any practical problem of interest.
Rather than compute the inverse of Z, an LU (lower-upper) factorization of Z can be performed, followed by a forward and back solve routine to solve for J. This requires less computational time than computing the inverse of Z and then computing J and is more accurate. The LU factorization of Z can be expressed as follows:ZJ=LUJ=Vwhere L and U are lower and upper triangle matrices. For example, for a 3×3 scalar matrix:
      [                                        Z            11                                                Z            12                                                Z            13                                                            Z            21                                                Z            22                                                Z            23                                                            Z            31                                                Z            32                                                Z            33                                ]    =            [                                    1                                0                                0                                                              l              21                                            1                                0                                                              l              31                                                          l              32                                            1                              ]        ⁡          [                                                  u              11                                                          u              12                                                          u              13                                                            0                                              u              22                                                          u              23                                                            0                                0                                              u              33                                          ]      Note that the elements of L and U are characteristic of the object and do not depend on the source function or excitation function. The L and U elements of an n×n scalar matrix are computed as follows:
            u      ij        =                  z        ij            -                        ∑                      p            =            1                                i            -            1                          ⁢                              l            ip                    ⁢                      u            pj                                          l      ij        =                  1                  u          jj                    ⁢              (                              z            ij                    -                                    ∑                              p                =                1                                            j                -                1                                      ⁢                                          l                ip                            ⁢                              u                pj                                                    )            Then, the J elements are computed by:
            x      i        =                  v        i            -                        ∑                      p            =            1                                i            -            1                          ⁢                              l            ip                    ⁢                                    x              p                        ⁡                          (                              forward                ⁢                                                                  ⁢                solve                            )                                                              with        ⁢                                  ⁢        i            =      1        ,    2    ,    …    ⁢                  ,          n      ⁢                          ⁢      and                  j      i        =                  1                  u          ii                    ⁢              (                              x            i                    -                                    ∑                              p                =                                  i                  +                  1                                            n                        ⁢                                          u                ip                            ⁢                              j                p                                                    )            ⁢              (                  back          ⁢                                          ⁢          solve                )                                with        ⁢                                  ⁢        i            =      n        ,          n      -      1        ,    …    ⁢                  ,    1  
In general, the system interaction matrix may be in block format for efficiency of organization and computation where each block corresponds to a group of unknowns. For example, each submatrix may express the interaction between each part (group of unknowns) of an object with itself or with other parts (group of unknowns) of the object. For the general case of an asymmetric block matrix, the LU factorization takes the form (for a 4×4 block matrix) where each block might correspond to tens to thousands of individual unknowns:
            [                                                                  [                Z                ]                            11                                                                          [                Z                ]                            12                                                                          [                Z                ]                            13                                                                          [                Z                ]                            14                                                                                          [                Z                ]                            21                                                                          [                Z                ]                            22                                                                          [                Z                ]                            23                                                                          [                Z                ]                            24                                                                                          [                Z                ]                            31                                                                          [                Z                ]                            32                                                                          [                Z                ]                            33                                                                          [                Z                ]                            34                                                                                          [                Z                ]                            41                                                                          [                Z                ]                            42                                                                          [                Z                ]                            43                                                                          [                Z                ]                            44                                          ]        ⁡          [                                                                  [                J                ]                            1                                                                                          [                J                ]                            2                                                                                          [                J                ]                            3                                                                                          [                J                ]                            4                                          ]        =          ⁢          ⁢                              [                                                    I                                            0                                            0                                            0                                                                                                          [                    L                    ]                                    21                                                            I                                            0                                            0                                                                                                          [                    L                    ]                                    31                                                                                                  [                    L                    ]                                    32                                                            I                                            0                                                                                                          [                    L                    ]                                    41                                                                                                  [                    L                    ]                                    42                                                                                                  [                    L                    ]                                    43                                                            I                                              ]                ⁡                  [                                                                                          [                    U                    ]                                    11                                                                                                  [                    U                    ]                                    12                                                                                                  [                    U                    ]                                    13                                                                                                  [                    U                    ]                                    14                                                                                    0                                                                                  [                    U                    ]                                    22                                                                                                  [                    U                    ]                                    23                                                                                                  [                    U                    ]                                    24                                                                                    0                                            0                                                                                  [                    U                    ]                                    33                                                                                                  [                    U                    ]                                    34                                                                                    0                                            0                                            0                                                                                  [                    U                    ]                                    44                                                              ]                    ⁡              [                                                                              [                  J                  ]                                1                                                                                                          [                  J                  ]                                2                                                                                                          [                  J                  ]                                3                                                                                                          [                  J                  ]                                4                                                    ]              =                  ⁢                  ⁢          [                                                                  [                V                ]                            1                                                                                          [                V                ]                            2                                                                                          [                V                ]                            3                                                                                          [                V                ]                            4                                          ]      where each Zij, is a matrix containing the terms of interaction between an ith source block (group of unknowns) and a jth observation block (group of unknowns). Note that the V and J may be matrices corresponding to multiple excitation vectors V and corresponding multiple source vectors J. For example, in monostatic scattering problems, one desires to determine J for each of a plurality of excitations, V. The factors L and U are determined as follows:
            U      ij        =                  Z        ij            -                        ∑                      p            =            1                                i            -            1                          ⁢                              L            ip                    ⁢                      U            pj                                          L      ij        =                  Dinv        jj            (                        Z          ij                -                              ∑                          p              =              1                                      j              -              1                                ⁢                                    L              ip                        ⁢                          U              pj                                          )      where Dinvjj=U−1jj. The forward solution for the blocks of Yi, with i=1, 2, . . . , n is:
      Y    i    =            V      i        -                  ∑                  j          =          1                          i          -          1                    ⁢                        L          ij                ⁢                  Y          j                    and the back solution for the blocks of Ji, with i=n, n−1, . . . , 1 is:
      J    i    =            Dinv      ii        (                  Y        i            -                        ∑                      j            =                          1              +              1                                n                ⁢                              U            ij                    ⁢                      J            j                                )  Once again, note that the D−1 operation can be computed equivalently by LU factorization, rather than by matrix inverse. That is, D can be LU factored, and then the left hand side can be forward/back solved for Ji.
For a symmetric block of Z submatrices, Z can be factored as follows:Z=U′TDU′.For example, consider:
      [                                        Z            11                                                Z            12                                                Z            13                                                            Z            12            T                                                Z            22                                                Z            23                                                            Z            13            T                                                Z            23            T                                                Z            33                                ]    =                    [                                            1                                      0                                      0                                                                          U                12                                  ′                  ⁢                                                                          ⁢                  T                                                                    1                                      0                                                                          U                13                                  ′                  ⁢                                                                          ⁢                  T                                                                                    U                23                                  ′                  ⁢                                                                          ⁢                  T                                                                    1                                      ]            ⁡              [                                                            D                11                                                    0                                      0                                                          0                                                      D                22                                                    0                                                          0                                      0                                                      D                33                                                    ]              ⁡          [                                    1                                              U              12                              ′                ⁢                                                                                                                          U              13                              ′                ⁢                                                                                                                            0                                1                                              U              23                              ′                ⁢                                                                                                                            0                                0                                1                              ]      where each Zij is a submatrix of Z. The diagonal entries of U′ are unity. The matrices Dii are the on-diagonal matrices shown below. Putting U′tDU′ into LU form, we first expand DU′ to obtain the form
      [                                        Z            11                                                Z            12                                                Z            13                                                            Z            12            T                                                Z            22                                                Z            23                                                            Z            13            T                                                Z            23            T                                                Z            33                                ]    =            [                                    1                                0                                0                                                              U              12                              ′                ⁢                                                                  ⁢                T                                                          1                                0                                                              U              13                              ′                ⁢                                                                  ⁢                T                                                                        U              23                              ′                ⁢                                                                  ⁢                T                                                          1                              ]        ⁡          [                                                  D              11                                                                          D                11                            ⁢                              U                12                ′                                                                                        D                11                            ⁢                              U                13                ′                                                                          0                                              D              22                                                                          D                22                            ⁢                              U                23                ′                                                                          0                                0                                              D              33                                          ]      Making the substitutionsUij=DiiU′ij Uii=Dii Lij=U′jiT=[Dii−1Uij]T=UjiTDii−1 we arrive at the standard LU form
      [                                        Z            11                                                Z            12                                                Z            13                                                            Z            12            T                                                Z            22                                                Z            23                                                            Z            13            T                                                Z            23            T                                                Z            33                                ]    =            [                                    1                                0                                0                                                              L              21                                            1                                0                                                              L              31                                                          L              31                                            1                              ]        ⁡          [                                                  U              11                                                          U              12                                                          U              13                                                            0                                              U              22                                                          U              23                                                            0                                0                                              U              33                                          ]      with factorization
      U    ij    =            Z      ij        -                  ∑                  p          =          1                          i          -          1                    ⁢                        U          ip          T                ⁢                  D          pp                      -            1                          ⁢                  U          pj                    and solve
                    X        i            =                        V          i                -                              ∑                          p              =              1                                      i              -              1                                ⁢                                    U              ip              T                        ⁢                          D              pp                              -                1                                      ⁢                          X              p                                            ;    Forward                      J        i            =                        D          ij                      -            1                          (                              X            i                    -                                    ∑                              p                =                                  i                  -                  1                                            1                        ⁢                                          U                ip                            ⁢                              J                p                                                    )              ;    Back  which completes the block symmetric matrix LU factor and solve expressions. Note that the inverse terms Dii−1 can be determined using LU factorization so that the inverse need not be computed directly.
Thus, LU factorization may be performed instead of a direct matrix inverse to solve the system matrix equation ZJ=V. However, although LU factorization of the system matrix equation is faster and more accurate than direct inversion of the system matrix, Z, determining the elements of Z and performing the LU factorization still takes a large amount of time (days) for objects that are electrically moderate or large in size and also requires large amounts of memory. Indeed, for many electrically large objects of practical interest, the solutions outlined above exceed the capacity of most if not all computing systems.
In recent years, another method has been developed for solving the system matrix equation ZJ=V. This method involves three basic steps: (1) spatially grouping the unknown sources, J, (2) computing a compressed form of the system interaction matrix, Z, and (3) performing iterative solutions to solve for J. This method avoids both matrix inversion and LU factorization, and further, compressing blocks of Z reduces the number of elements of each Z block that must be computed and stored. The method exploits the fact that when the unknowns are spatially grouped, then off-diagonal sub-matrix blocks of Z are sparse. This means that the interaction physics contained in a Z matrix block describing the interaction of a group of, for example, 3000 unknowns with a distant group of 3000 unknowns is a block matrix of 3000×3000 (9 million complex numbers) that can be expressed with far fewer numbers. Mathematically, this 3000×3000 block is of low rank and can be compressed with significantly less memory storage. Moreover, since the currents are solved iteratively, a full matrix inversion or LU factorization need not be performed.
Spatial grouping of unknowns involves ordering unknowns according to distance between a point selected as a reference point and other points on the object. This results in a system matrix with low rank block submatrices. The blocks furthest from the diagonal of the system matrix are generally the blocks of lowest rank. Because the blocks are of low rank, they may be partially computed in a compressed form. One method of computing a matrix Z in compressed form is described in: Mario Bebendorf, “Approximation of boundary element matrices,” Numer. Math. (2000) 86:565:589, which is incorporated herein by reference to “Bebendorf” for all purposes. This method, known as the Adaptive Cross Approximation (ACA) technique, is also described in: Stefan Kurz, “The Adaptive Cross-Approximation Technique for the 3-D Boundary-Element Method,” IEEE Transactions on Magnetics, Vol. 38, No. 2, March 2002, which is also incorporated herein by reference to “Kurz” for all purposes.
The ACA technique may be outlined as follows. Consider a matrix A to be approximated as the outer product of Au times Av where Au is a column dominant matrix and Av is a row dominant matrix. Any row or column of A can be computed from a subroutine, GetRowOrColumn. Thus, an arbitrary row of A is first computed. This row is a pivot row, kRowPivot, v0. Then select the column of this row containing the largest element. Call this column kColPivot. Divide the row, v0, by this largest element, so that all the elements of row v0 are now unity or less. Now compute the kColPivot u0 column from the subroutine, GetRowOrColumn. The outer product, u0v0, is the first approximation to the matrix, A, A0=u0v0. Let an error matrix be E0=A−A0. The next row and column approximation, u1, v1 is computed, just as described for u0, v0, from the error matrix, E0, and the next approximation is A1=u0v0+u1v1. The next error matrix is E1=A−A1. Thus, the process is continued until the error matrix is sufficiently small. If k terms are needed for convergence, the k uk column vectors comprise the matrix Au and the k row vectors vk comprise the Av matrix. If the matrix A is low rank, k will be much less than m or n original dimensions of matrix A.
Thus, using ACA, a compressed form A can be computed as:[Ac]=[Au][Av]where Au and Av are column and row dominant matrices, respectively of size (m×k) and (k×n). The vectors u and v comprising Au and Av are computed recursively using the method described in Bebendorf as shown in Kurz. Bebendorf has shown that to compute this approximation to A, only k rows and columns of A are needed. If A can be approximated with k << (m or n), then a significant reduction in matrix fill time can be accomplished. Thus, compression results from the fact that only a relatively small number of rows and columns of A need to be computed to obtain the resultant compressed matrix, Ac. A measurement of the error can be obtained as:
                                                                    A              -                                                A                  c                                ⁡                                  (                  k                  )                                                                                                  A                                      =        ɛ                            (        2        )            If matrix Ac=Au Av represents A to this tolerance, we say that A has rank k for tolerance ε. Computing only a relatively small number of rows and columns of A, this tolerance error can be on the order of 10−4, for example.
In this method, when the compressed form of a Z block is computed as ZuZv, the solution for J is computed iteratively. In the iterative process, a first approximation, J1, is substituted into the matrix equation and the resultant excitation, V1, is computed. An error vector e=V−V1 is computed. This error is used to compute a second, hopefully more accurate, approximation, J2, and the process is repeated. This continues until the error becomes sufficiently small. Methods for implementing iterative solutions are well known in the art.
Several problems exist with the compression and iterative solve method just described. First, although the method of computing a compressed system interaction matrix is substantially faster than computation of the full system interaction matrix, computation of the iterative solution still takes a very long time for problems of practical interest. Second, for many problems of practical interest, the iterative solution converges very, very slowly or not at all, such as for cavities which are of great interest in stealth design problems. And for each RHS (right hand side) vector V a new iterative process must be done. For backscatter problems for electrically large bodies with many RHS to describe the pattern, iterative approaches require an inordinate number of RHS.
Thus, all of the above-described methods for solving a system matrix equation have the limitation that computational time is excessively long and required memory storage is excessively large. What is needed is a method that substantially reduces computation time and memory storage while producing accurate results.